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Abstract 

In this paper the optimal transport and the metamorphosis perspectives are combined. For a pair of given input images 
geodesic paths in the space of images are defined as minimizers of a resulting path energy. To this end, the underlying 
Riemannian metric measures the rate of transport cost and the rate of viscous dissipation. Furthermore, the model is 
capable to deal with strongly varying image contrast and explicitly allows for sources and sinks in the transport equations 
which are incorporated in the metric related to the metamorphosis approach by Trouve and Younes. In the non-viscous 
case with source term existence of geodesic paths is proven in the space of measures. The proposed model is explored on 
the range from merely optimal transport to strongly dissipative dynamics. For this model a robust and effective variational 
time discretization of geodesic paths is proposed. This requires to minimize a discrete path energy consisting of a sum 
of consecutive image matching functionals. These functionals are defined on corresponding pairs of intensity functions 
and on associated pairwise matching deformations. Existence of time discrete geodesics is demonstrated. Furthermore, 
a finite element implementation is proposed and applied to instructive test cases and to real images. In the non-viscous 
case this is compared to the algorithm proposed by Benamou and Brenier including a discretization of the source term. 
Finally, the model is generalized to define discrete weighted barycentres with applications to textures and objects. 


Introduction 

In the past two decades concepts from finite dimensional classical geometry have been successfully transferred to infinite¬ 
dimensional spaces, where shapes are contour curves of geometric objects, surfaces, image intensity maps, or probability 
densities. These concepts have a continuously increasing impact on the development of novel computational tools in 
computer vision and imaging, ranging from shape morphing and modeling BKMP07I . and shape statistics, e.g. | |FLPJ04 |, 
to texture analysis I RPDB12 1 and computational anatomy IIBMTY021 . Three particularly influential approaches on the 
space of image maps are linked to optimal transportation IIMonS 11 iBBOO. ' Z YHT071 [ PPO14 1, the flow of diffeomorphism 
BAm661 lDGM98al and metamorphosis [MY01 TY05|. Here, we combine these three approaches and explore properties 
of the resulting image manifold. Thus, in what follows we briefly review the underlying concepts. 

Optimal transport and its application in imaging. The problem of optimal transport is introduced in the seminal 
work of Monge I IMon81 1 in 1781. In I Kan421 |Kan481 Kantorovich proposes a relaxed formulation of Monge’s problem 
which gives rise to the Wasserstein distance considered in this paper. Let ( D. d) constitute a metric space. The 2- 
Wasserstein distance between two probability measures ^a, dB £ P(-D) is defined by 

W(ha,^b) 2 ■= min d{x,y) 2 dir(x,y). (1) 

7rer Jdxd 

Here T (ha, dn) denotes the set of all probability measures on D x D with marginals ha and /is with respect to x and 
y , respectively. For an introduction to optimal transport and the Wasserstein distance we refer the reader to the reviews 
IIAmb03l IEva99l IViKGl IAGS06I I Vil08l , 

Being a distance function applicable to very general measures (continuous and discrete measures) the Wasserstein 
distance has an increasing impact on robust distance measures in imaging iRTGOOt BS051IPFR121 [BFS121 . As the 
Wasserstein metric is defined for arbitrary probability measures with finite second moment it allows to measure distances 
between absolutely continuous measures with respect to the Lebesgue measure as well as concentrated measures. With 
the increase of the complexity of applications efficient numerical computation of ([3]) became increasingly important. In 
that respect Benamou and Brenier propose an alternative formulation of the quadratic Wasserstein distance using the 
perspective of the underlying flow of a density 6 with Eulerian velocity v and expressing the transport in terms of a 
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constrained flow BBB001 . That is, one asks for a minimizer of the path energy 


£[0,v\ 


e\v\ 2 dt 


D 


( 2 ) 


for a density function 0 : [0,1] x D —» K. and a velocity field v : [0,1] x D —>• subject to the transport equation 

d t 0 + div( 0 i>) = 0 and the constraints 0 ( 0 ) = 9 a and 0 ( 1 ) = 0 B . Then the minimal energy is indeed the squared 
Wasserstein distance between measures ha and //« with corresponding densities 9 a and 9 B respectively. Their algo¬ 
rithm has been immensely influential in the numerical computation of the Wasserstein distance and gradient flows related 
to it, see e.g. 1BCW101 IDMMlOl IPPQ141 . In IPPQ141 , for instance, a proximal point algorithm for the solution of 
Benamou-Brenier’s formulation (J2| is derived and applied for the computation of Wasserstein geodesics between two 
image densities. Alternatively, if D C M. d is a strictly convex domain, d is the Euclidean distance, and /i.a and fin are 
absolutely continuous measures with densities 9a and 9 B , respectively, then one has 


>V(/m,M.b ) 2 = min / \(f>(x) - x\ 2 9 A {x)dx, (3) 

where denotes the push forward of the measure // 4 under the mapping 0 and for a diffeomorphism 0 the constraint 

can be expressed as (det Dcf>)9 B o <fi = 9a- In I AHT03 I an initial mass preserving transport map li Mos65 1 is created 
and then an explicit time stepping scheme is employed to compute the optimal map from a modified formulation of <[3j 
where the constraint is linearized. In [ HRTlOl l the authors pick up formulation |3]i as well and use a sequential quadratic 
programming method for its optimisation. Moreover, in I CWVB09 1 the authors propose a gradient descent for the dual 
formulation of ([3} and show its use for image registration and warping. In K LR05I ISAKlOl a damped Newton method 
is used to compute a solution rp of the Monge-Ampere equation (which is the equality constraint on 0 = VTO and sub¬ 
sequently the optimal transport map 0. Finally, let us mention that in 1 ISS 13all the authors propose another interesting 
numerical algorithm for the efficient computation of the Wasserstein distance that is based on an extension of the auction 
algorithm. The latter optimizes the cost functional in the Wasserstein distance only on a sparse subset of possible assign¬ 
ment pairs still guaranteeing global optimality. Various alternative computational approaches for the Wasserstein distance 
exist, e.g. HDG06I IQbeOSl IPPO141 . 

In terms of imaging application the Wasserstein distance has been employed in the context of image and shape classi¬ 
fication, segmentation, registration and warping, and image smoothing. In [ RTG00 I the Wasserstein distance is used as a 
distance on the images directly, interpreting images as discrete measures. In other approaches the Wasserstein distance is 
used on image histograms or points clouds (which could be feature vectors of images), see e.g. BGD04IILQ071 . In the con¬ 
text of image segmentation and classification similar approaches are used, see, e.g. [CEN07 NBCE09l lPFRT2llQJBS12l . 
In I RPDB1 2 I a discrete Wasserstein distance is employed for the computation of barycentres of discrete probability distri¬ 
butions with applications to texture synthesis and mixing. Thereby the original Wasserstein metric is replaced by a sliced 
version over one-dimensional distributions. In I RPC 10 1 the same approach is used on clouds of geodesic shape descriptors 
as a similarity measure to discriminate between different shapes in 2D and 3D shape retrieval. The Wasserstein distance is 
also used in the context of contrast and colour modification in, e.g. HRP11IIFPR+I31 . In BZHT03IIHZTA041 the quadratic 
Wasserstein distance is considered to define a rigorous distance between images, applied to non-rigid image registration 
and warping. As a distance function for shapes the Gromov-Wasserstein distance is introduced in a series of works by 
Memoli IIMem071 IMem 111 . Here, shapes are modelled as compact metric spaces and the Gromov-Wasserstein distance 
is computed on isometry classes of each space. The Gromov-Wasserstein differs from the Wasserstein distance G as it 
assigns a cost to pairs of transport assignments. In the context of surface dissimilarity measurement the Wasserstein dis¬ 
tance is applied to metric densities on the hyperbolic disc representing conformal mappings of different surfaces l lLDl 1 1. 
In 1ISS13bl the authors create a convex shape-prior from a modified Gromov-Wasserstein distance. Their approach can 
be used for image segmentation problems in which prior shape knowledge on the objects that should be segmented can 
be provided in terms of a template shape. This approach is modified in 1SS 13cll using the quadratic Wasserstein distance 
as a regularizer between learned reference shapes and the segmentation. In I BFS12 I the Wasserstein distance is used 
as a data fidelity term in a generic regularization approach and applied for image density estimation and cartoon-texture 
decomposition. Moreover, a review of the use of geodesic methods and in particular optimal transportation in computer 
vision can be found in I PPKCIO I. 

Flow of diffeomorphism. The physical modeling of viscous flow involves dissipation as an integrated measure of 
local friction. Arnold I Arn66 ) AK98 1 proposes to study viscous flows from the perspective of a family of diffeomorphisms 
(< 5 f>(f)) te [o,i] : D —> D which describe the transport of densities, e.g. image intensities, along particle paths (<(>(i, a;))te[o,i] 
for x £ D. This concept is picked up in vision by Grenander and coworkers BGre8111DGM98b l. As a Riemannian metric 
Q dlS [-, ■] one considers the rate of viscous dissipation induced by the Eulerian flow velocity v(t) = (j){t) o ((> - 1 (t) in a 


2 






















































multipolar fluid model (cf. Necas and Silhavy HNv91iO . Here, the Eulerian motion field v is considered as a tangent vector 
on the manifold of diffeomorphisms. The resulting Riemannian metric one obtains G‘ uff [v. u] := f D L[v, v\ dx, where 
L[v, v ] = Ce\v\ : e[u] with e[v\ = (Vw + (Vx) T )/2, and a deduced path energy £ dlff [</>] = f* [u, u] d t as an action 
functional on flows encodes the total accumulated dissipation on the domain D and on the time interval [0,1], To study the 
warping of two image intensity functions 9 a, 6 b '■ D —► R. for which there exists a diffeomorphism <fi with 9b = 6a ° 4> 
a flow minimizing the energy £ dlff subject to the constraints <j>(0) = II and 6b = 6a ° 4>{ 1 ) defines a geodesic path 
]i with 9{t) = 6a o <j> 1 (t) in the spaces of images connecting 6a and 9b- If we aim at deriving a Riemannian 
distance directly between images via the flow of diffeomorphism approach, then a motion field v can be viewed as a 
representation of an image variation. Obviously, different motion fields might represent the same image variation. Hence, 
the corresponding equivalence class v is considered as a tangent vector on the image manifold and the associate metric is 
now given by 

f/ dlff [u, v) ■.= min f L[v,v]dx. (4) 

v ^ v Jd 

Consequently, the path energy on a path (0(i)) tg [ 0 i i] reads as 

£ diff [0] = f G diS [v,v]dt, (5) 

Jo 

which one minimizes over all image paths with 0(0) = 6a , 0(1) = 0s- Here, we assume that there is at least one 
path of finite path energy connecting 6a and 9b- In medical applications IBMTY02 I each diffeomorphisms represents a 
particular anatomic configuration of an anatomic reference structures. For more details we refer to HDGM98al 1BMTY051 
iJMOOl ImTY021 . 

Metamorphosis. A one-to-one correspondence of image grey values in warping applications is frequently not re¬ 
alistic. The metamorphosis approach offers a suitable generalization of the flow of diffeomorphism concept. It is first 
presented by Miller and Younes [MYOlj. A rigorous analytical treatment is due to Trouve and Younes I TY05 I. In ad¬ 
dition to the transport of image intensities along motion paths the variation of an intensity value along a motion path is 
allowed and reflected by an additional term in the energy. This term measures the integrated squared material derivative 
z = Of f) + V0 ■ v. From a geometric perspective a pair of material derivative z and motion velocity v represent a variation 
of an image 0. Hence, an equivalence class (z, v) of all pairs which generate the same image variation is considered as 
a tangent vector on the image manifold. A Riemannian metric Q meta [(z, v), (z, u)] acts on these tangent vectors and the 
associated path energy along an image path ( 0 (£))te[o,i] is given by 

£ meta [0] = f 1 g meta [M, (z^)} dt. (6) 

Jo 

As an example for the underlying Riemannian metric we obtain 

£ meta [M, {^v)] = min f L[v, v] + jZ 2 dx , (7) 

(z,v)e(z,v) Jd ° 

where the first three terms in the integrant retrieve the metric from the flow of diffeomorphism approach and encode the 
induced viscous dissipation, whereas the last term penalizes temporal changes of intensities along motion paths. 

In this paper, we combine the optimal transportation approach with the metamorphosis approach. Thereby, in addition 
to the transportation cost we take into account a density variation of the transported measure and viscous dissipation. The 
paper is organized as follows: First we present our generalized image transport model in Section[T] In Section[2]we prove 
existence of geodesics in the non-viscous case. Then we propose a variational time discretization of the full model in 
Section[3] prove existence of time discrete geodesics in Section[4]and describe a fully discrete solution scheme in Section 
[5] Furthermore, we consider in Section[ 6 ]the algorithm used inl lBBOOl to compute for comparison reasons geodesics in 
the purely non-viscous case. Finally, we generalize in Section[7]our model to discrete weighted barycentres and apply it 
to textures and objects. 

1 The generalized image transport model 

In this section we will discuss the generalization of the optimal transport model in image warping and blending. These 
generalizations are motivated by two observations in applications: 

- Frequently, objects or structures in images, which are in correspondence and are expected to be matched via the transport. 
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have different masses. From a global perspective the assumptions that images are considered as probability distributions 
is too restrictive. Indeed, the latter requires in advance contrast modulation, which is somewhat artificial. In the classical 
optimal transport model local mass differences lead to artifacts, where a local mass surplus has to be deposited elsewhere 
without any structural correspondence. We will no longer enforce the source free transport equation and explicitly incor¬ 
porate a source term in the path energy which measures density modulation. 

- Different from the flow of diffeomorphism approach the optimal transport maps are not necessarily homeomorphisms. 
On the other hand in many applications one is interested in topological consistency and at the same time the physical 
background of the application might suggest to incorporate a dissipative term in the path energy. Hence, we combine the 
classical transport cost model with a weighted viscous dissipation model. 

To this end we first recall the formulation ([2) of the Wasserstein distance proposed by Benamou and Brenier l iBBOO l. 
In what follows we restrict to a bounded domain D C (d = 2, 3) with Lipschitz boundary. Now, we allow for a source 
term z : [0,1] x D — > 1 in the transport equation defined for given image intensity 0 : [0,1] x fl — y R and transport field 
v : [0,1] X D —> R d as 

z := d t 9 + div( 0 u). ( 8 ) 

Furthermore, we pick up the model for the viscous dissipation in |5]i and obtain as a new path energy 

£s n [9]= [ Gs^[(z,v),(z,v)}dt. (9) 

Jo 

with 

g 5 ^[(z,v),(z,v)\= min_ f 9\v\ 2 + \z 2 + <yL[v, v] da:, (10) 

{z,v)G(z,v) JD ° 

which we minimize subject to <[ 8 ]> and the constraints 0(0) = 9a and 0(1) = 9b- Here, 9a and 9b are the given 
input images and the (z, v) is the equivalence class of pairs of a source term and a transport field, which are consistent 
with the transport equation ([S} for given image intensity 0. The involved local rate of viscous dissipation is given by 
L[v, u] = |(tre[u ]) 2 + ^tr(e[u] 2 ) + e\D m v\ 2 , where e[u] = |(Vv + Vu T ), m > 1 + f and A, fi, r] > 0. (The first two 
terms represent the viscous dissipation of a Newtonian fluid and the higher order terms reflect a multipolar viscosity). As 
in I TY05 1 and similar to 1AGS06 1 the condition 9 t 0 + div(0u) = z has to be understood in weak form 


rjz dx di = — 


o Jd 


(d t ri + v ■ \7r/)0 dx d t 


o JD 


for all rj £ C™{{0,1) x D). 

The first term in the metric is the classical transport cost rate, the second terms reflects the source term which measures 
the density modulation of the image intensity and the last term is the dissipation rate based on a multipolar viscous fluid 
model. Let us emphasize that for general non divergence free motion fields z does not coincide with the material derivative 
as in the metamorphosis model IITY05I . We suppose that 7 > 0 measuring the impact of viscosity and 0 > 0 is a penalty 
parameter weighting the impact of density modulation on the metric and the path energy. In the formal limit S —> 00 
and for 7 = 0 we retrieve the standard transport cost. Given the path energy, we can define a Riemannian (generalized 
Wasserstein) distance Ws tl \9A-> 0b] of two images 9a and 9b as 

WsA 9 a,0b] 2 = min £ a , 7 [0]. (11) 

(0(*))t£[ 0,1] 

0(O) = 0 A , B(1) = 6 b 

2 Existence of geodesics for the non-viscous model (7 = 0 ) 

In this section we study existence of minimizers 0 of (|9]> in the non-viscous case, that is for 7 = 0 . In order to give 
a rigorous proof, it will be necessary to reformulate the formal problem (|9| as a problem for measures rather than for 
densities. In particular, it will be crucial to treat the singular parts of the measures in an appropriate way. 

Following IBB00I , it will be useful to replace the velocity variable v by the momentum variable w = 9v. Therefore, 
the function 0|x| 2 appearing in the path energy ( |T0| will be replaced by uj 2 /0. The joint convexity of this function will 
play a crucial role in the sequel. 

The argument presented here is a modification of the argument in I DNS09 1 and our presentation follows this latter 
work very closely. Some additional arguments are needed to deal with the possibly varying total mass. On the other hand, 
some simplifications can be made, since we work on a bounded spatial domain instead of the whole space R d . 
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First we reformulate the action functional ([TO]). Let Q he a bounded domain in Euclidean space and fix a reference 
measure Jz? £ ■// + (f-l). In our application, the domain Cl will either be the spatial domain D or the space-time domain 
[0,1] x D, and Jz? will be the corresponding Lebesgue measure. 

Let p £ ^f/ + (Cl), v £ (Cl' 1 M. d ), and ( £ l). The Lebesgue decomposition of these measures with respect to 

Jz? is given by 

H = 0Jz? + , v = u>Jz? + u 1 - , ( = zd£ + £-*- . 

Let now .16' £ ./// + (f.l) be such that // L . v 1 and are absolutely continuous with respect to Jz? 1 (take for instance 

^ + \v x \ + |(- L |). Then we may write 

/i x = e ± £ ,± , ^ = w ± £ ,± , c 1 - = z ± ^ ,± . 

As in the Benamou-Brenier formulation of the 2-Wasserstein distance we consider the function f : [0, oo) x K d —> [0, oo] 
defined by 

f 0 0 = 0 and w = 0 , 

0 > 0 , 

[ +oo 0 = 0 and (u/ 0 . 


Note that <\> is lower-semicontinuous, convex and 1-homogeneous. The action functional V : x x 

( 12 ) —* [ 0 , +oo] that we are interested in is given by 


V, C) := X>bb(/L v) + \ V z(C) , 


T^bb (AT v) ; — 

%(C) := 


(f>(9, w) dJz? 


4>(0 J 


,w 


)dJz* J 


f In ~ 2 d -^ C x = 0 , 

\ +oo C X / 0 ■ 


Since <j> is jointly 1-homogeneous, the definition of 2 ?bb does not depend on the choice of Jz?- 1 . The same is true for 
T>z, since we may write T>z{ C) = J n z 2 dJz? + J Q ip(z ± )d^f- L , where if) : K. —► [0, +oo] is the 1 -homogeneous function 
defined by tp(0) = 0 and f>(r) = +00 for r / 0. Sometimes it will be useful to write V n instead of D in order to 
emphasize the domain 12 . 

The following result is an immediate consequence of general lower-semicontinuity results for integral functionals on 
measures 1 AB881 AFP00 I. 


Proposition 2.1 (Lower semicontinuity of the functional V). Consider weak*-convergent sequences of measures 


At n * /it € JZ + {Ci) , v n -^* v £ Jf{Cl\ R d ) , C-x^(12) . 


Then we have V(n, v, () < liminf^oo T>(fj, n , u n , ( n ). 

Proof. The result follows, since both 2 ?bb and T>z satisfy the assumptions of [ DNS091 Theorem 2.1], O 


The following crucial lemma is a special case of IDNS091 Proposition 3.6], 

Lemma 2.2 (Integrability estimate). Let a t £ ^ + (f l) and v £ ,/CT, (12; R d ). For any Borel function r/ : 12 —> K + we have 


r](x)d\a\(x) < (v BB (n,u)Y f f r) 2 d^t 


Proof Set S := {x £ Cl : rj(x) > 0}. Using the scalar inequality sfab + \fcd < y/a + cy/b + d which holds for 
a, b,c,d> 0, we obtain 


rj(x)d\v\(x) = f p|u;| dJz? + f r]\'w^\dd£ 1 ^ 
Js Js 


< ( / <j>(9,w) dJz? ) ( / ry 2 0dJz? 




f(9 ± ,w ± ) dJz? J 


ti 2 9 ± d^? x 


< ( / 0(0, w) dJz? + / cj)(9 ± ,w ± )dJf J 
s Js 


if9d^f+ I tf9 ± di? x 


< 


(v m (ax, 1 /)) 2 (^ 2 dd 


□ 
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Let us now introduce the modified continuity equation. 


Definition 2.3 (A continuity equation without conservation of mass). Let t H>■ p t be weak*-continuous in (I)), let 
t i—>- be Borel measurable in „•#( D;WL d ), and let t t-A ft be Borel measurable in ^fT(D). We say that the triple 
(p t , n t , Ct)te[o,il satisfies the continuity equation (and write (p t , v tl Ct)te[o,i] £ C£[0,1 ]) if 

1. the following integrability conditions hold: 

[ Wt\(D)dt < oo , [ |C t |(D)df < oo ; 

Jo Jo 

2. the modified continuity equation dtp + div(i/) = f holds in the sense of distributions, i.e., for all space-time test 
functions rj £ Cg((0,1) x D) we have 

f f d t p(t,x)dp t (x)dt + f f Vp(t,x) ■ du t (x)dt+ f f p(t, x) df t {x) dt = 0 . 

Jo Jd Jo Jd Jo Jd 

A standard approximation argument (see II DNS091 Lemma 4.1]) shows that solutions to C£[ 0,1] satisfy, for all 0 < 

t 0 <h< 1, 


ID 


v(h,x)dp tl (x) ~ / ??(to ,x)dp to (x) = 


ID 


’ to J D 


d t p(t,x)dp t (x)dt - 


/to J D 


V?7(f, x) ■ dv t {x)dt 


( 12 ) 


/to J D 


ri(t,x)d( t (x)dt 


for all space-time test functions p £ C 1 ([0,1] x D). In particular, taking p(t, x) = 1, it follows that the increase of mass 
is given by 


PtA D ) ~ Bto(D) = / (t(D)dt . 


(13) 


J to 

We are now in a position to rigorously define the extended distance Ws that was formally introduced in {13- 

Definition 2.4. For pa,Pb £ + (D) we define Ws(pa, Pb) £ [0, +oo] by 

Ws(pa,Pb) ■= j (^J T>(p t ,v t ,( t )d?J : (p t , u t , Ct)te[o,i] £ C£ [0,1] , p 0 = Pa , Pi = Pb^ ■ (14) 

The following theorem is the main result of this section. 

Theorem 2.5 (Existence of geodesics). Let S £ (0, oo) and take pa, Pb £ ^ f/' + (D ) with Ws(pa, Pb) < oo. Then 
there exists a minimizer (p t ,u t , Ct)te[o,i] that realizes the infimum in ( | 14| >. Moreover, the associated curve (/Z t ) te [o.i] is a 
constant speed geodesic for W$, i.e.. 


Ws(p s , Pt) = |s - t\Ws(pA, Pb) 

for all s,t £ [0,1]. Furthermore, we have the alternative characterisation 

Ws(pa, Pb) ■= jrd j J s/V(p t ,u t ,Ct)dt : (pt,Vt,Ct)te[ o,i] £ C£[ 0,1] , po = PA , Pi = Pb j • 

Proof The existence of a minimizer is an immediate consequence of Proposition |2.6| below. The remaining statements 
follow by standard arguments, see I 1DNS091 Theorem 5.4] for details. □ 


Let us now state and prove the main ingredient for the proof of Theorem 2.5 We write f 0 S t 0 pt dt to denote the 
measure p on [0,1] x D satisfying 


' [ 0 ,l]x.D 


p(t,x)dp(t,x) = / / p(t,x)dp t (x)dt 

Jo J D 


for all 77 £ C([0,1] x D). 
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Proposition 2.6 (Compactness for solutions to the continuity equation with bounded action). Suppose that 

(m".C") te(o,i) G Cf[0 ,1] satisfy 

(Al) Mi := svtp n p,Q(D) < oo; 

(A2) M 2 := sup„ fg X>(/#, i/ t n , Q) df < oo. 

Sef v n := fg 5t®ft dt £ .#([0,1] x D\ K d ) and Cf := fg St®(pdt £ --#([0,1] x I?). Then, there exists a subsequence 
(again indexed by n) and a triple ( p t , r't, Ct)te(o,i) G C£[0, 1] such that 

1. pf —*■* pt in (D) for all t £ [0,1]; 

2. v n * v in y/([0, 1] x D; R d ); 

3. C" -"* C in -#([0,1] x D). 

Moreover, for the above subsequence 


f T>(p t ,u t ,Q)dt < liminf f T>{p?, uf, Q) df . 
Jo "-*- 00 Jo 


Proof. In view of (A2), we first observe that 

M 3 := sup f |£"|(£))df = sup f f \z^(x)\dxdt 
n Jo n Jo JD 

Therefore, ( p~3] > yields the uniform bound 

p n t(D) < Po(D) + f |CIP)ds < Mi + M 3 
Jo 

for all n. Moreover, Lemma |Z2] implies that 

Wt\(D) < Vtf(D)V BB (p?,v?) , 


< oo . 


(15) 


hence by the Holder inequality we obtain 


r K\(D) 2 dt < [ p?(D)V(p?,v?,Q)dt < M 2 (Mi +M 3 ) , 

Jo Jo 

which shows that the maps {t H>• \vfj \(£))}„ are uniformly bounded in L 2 ( 0 , 1 ), hence uniformly integrable. 

Since |^ n |([0, l]xfl)< (fg \v’f’\(D) 2 dt)^, the measures {a n }„ £ .M ([0,1] x D;M d ) have uniformly bounded total 
variation on [0,1] x D, hence we can extract a subsequence that converges weakly* to some measure v £ . M ([0.1] x 
D; R d ). The uniform integrability of (f i—>■ \vf\(D)} n implies that the image measure of u under the mapping (t, x ) i —> t 
is absolutely continuous with respect to the Lebesgue measure on [0,1]. Therefore, the disintegration theorem (see, e.g., 
I AGS06] Theorem 5.3.1]) allows us to write v — f n S t ® v t df for some family of measures {^t}tg[o,i] G ( D ; K d ). 

Fix 0 < r < 1, take q £ C^(D), and set £(f, x) := S/q(x)x\o, T ] (i)- Although £ is discontinuous, general approxima¬ 
tion results (see I AGS06 Proposition 5.1.10]) imply that 


Vq(x)dvf(x)dt = 


id 


' [0,1] xD 


£(t,x)dv n (t,x) 


' [0,1] X D 


^(t,x)du(t,x) = 


S/q(x)dv t (x)dt . (16) 


ID 


Let us now consider the term involving (f, which is treated similarly. For all n and a.e. t £ [0,1] we use (A2) to 
conclude that Cf = ■ Therefore we obtain 


\<f?\(D) 2 dt 




1 (x) | da; 


2 

df < oo . 


As above, we infer that the mappings {f i—>■ |£"| (D)} n are uniformly integrable, and that there exists a subsequence of 
{£"}„ that convergence weakly* to some measure ([0,1] x D ). By the disintegration theorem we may write 
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C = f 0 S t CS> (t d t for a family of measures {Ct}te[o,i] £ Set £(t,x) := t]{x)x[o,t} (£)• Arguing as above, we 

obtain 


/ T](x)dCt(x)d t= Z{t 1 x)dC{t,x) 

Jo Jd J[o,i]xd 


/ |(f,x)dC(f,x) = / / r/(x)dCt(x)df . (17) 

/ [0,1] x d 7o 7 _d 


We are now in a position to obtain subsequential convergence of {//,J' }„. Indeed, it follows from ([T2|> that 


/ v(xWt(x) = / r](x)d$(x) + / / Vtj(x) • d^ n (a:)df + / / r?(a;)dCrO)d£ 

Jd Jd Jo Jd Jo Jd 


Moreover, (Al) implies that there exists a measure y ,o € ^/£ + {D) such that ylfi —/io (after passing to a subsequence). 
In view of © and ( fUf the latter equation implies weak*-convergence of {/x"} n to some measure f,i T for every r £ [0,1]. 
It is readily checked that (y t , v t , Q) t £ C£[ 0,1]. 

It remains to prove ( fT5j ). For this purpose we write y, = f* St (g> yt d t and y n = St <E> p." df. It is straightforward to 
check that y n converges weakly* to y, in ^+([0,1] x D). Now the result follows by observing that 

[ V D (fi t ,v t Xt)dt = V^x D ^,v,0 , 

Jo 

and applying Proposition |2. 1 |to pi 0 ’ 1 ! x D _ 0 


3 A variational time discretization 

In what follows, we derive a time discrete approximation of the energy ([9]) and thereby a variational approach for the 
definition of geodesic paths. We refer to I WBRS11 1 for the general concept and to I RW14 I for the numerical analysis in 
the context of shape spaces which are Hilbert manifolds and in I BER15 I a variational time discretization of geodesics in 
the metamorphosis model is discussed. 

As a motivation let us briefly present a toy model in finite dimensions. On a smooth m-dimensional manifold M 
embedded in (to < d) we consider the simple energy F[y, y\ = \y — y\ 2 which reflects the stored elastic energy in 
a spring spanned between points y and y through the ambient space of M in IX f/ . The smoothness of M implies that 
F[y, y] = distal (y, y) 2 + 0(dist^vi (y, y) 3 ), where dist^ (y, y) denotes the Riemannian distance between y and y. Hence, 
we can approximate the path length of a smooth path (y(£))te[o,i] y i a sampling y\. = y{ and then evaluate the discrete 
path energy 

K 

E A [y 0 , ...,y K ]=K^2 \y k - y k _ i| 2 , 

k= 1 

such that E a [j/o, • ■ • ,Vk\ converges to f [(y(£))te[o.i]] = fo ll /(^)| 2 dt for K —> oo. Here, we use that is an 

approximation of the velocity y(k/K) where r = A is the time step size of our discretization on the time interval [0,1]. 
In fact, based on the approximation F of the squared distance distjvt, which is easy to implement, we obtain an effective 
approximation of the Riemannian path energy £. Correspondingly, we call a minimizer (j/q, • ■ • ,Vk) of the discrete 
path energy for fixed yo and y k a discrete geodesic. We refer to l iRWl4 „ IBER15 I for a detailed discussion, why the 
discrete path energy instead of the discrete path length is the right concept to compute discrete geodesics. In particular, 
I'-convergence of the discrete path energy is proven in case of the metamorphosis model in I BER15I and under suitable 
assumptions in the context of Hilbert manifolds in IIRW141 . 

Now, we ask for a similar time discrete approximation of the continuous path energy Es rj defined in (|9|. To this 
end, we consider a discrete path {0$, ,9k) in the space of image intensities with 9 k £ I for k = l,... ,K with 
I := /j 2 (D. K> () ) and ask for a matching functional F on consecutive pairs 9, 9 of image intensities. In fact, this 
matching functional should reflect time discrete counterparts of all three ingredients of the metric Q$ n and the induced 
continuous path energy namely the transport cost, the viscous dissipation and the source term. Like in the original 
Monge problem we take into account deformations 0 in a suitable space A of admissible deformations, to be defined later, 
and optimize for given 9,9 a suitable functional F D [9, 9. ■] over all admissible deformations to define the value of the 
matching functional F[0, 9\, i.e. 

F [0,0\ = inf F D [9,9,(/)}. 
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With the matching functional at hand we then define the discrete path energy summing over applications of the matching 
functional to consecutive pairs of image intensities Ok- 1 and Ok of a discrete path (Oq, ... ,0k) and get 

K 

E?„[0o,...,0 K ]=K'£ i F[9 k - 1 ,0 k \. (18) 

k =1 


Thus, the resulting time discrete approximation of the squared Riemannian distance is given by 


W^[0a,0 b Y= n min_ Ef„[9 0 ,. 


So 

Oq =6 A: 


ex 

K—0 


'K\ • 


(19) 


Here, we assume 0 a , 0 b € X. In what follows we list now the appropriate components of F reflecting the different 
ingredients of the continuous path energy. 

Approximation of the transport cost. To approximate the first term in the metric ([TO]) we make use of the equivalence 
of the original Monge problem and the Benamou Brenier formulation I BB00 I of optimal transport and define 

Fl p JM= / \^-A\ 2 0dx. (20) 

J D 

Here, A-Jl is an approximation of the transport velocity with II being the identity deformation. 

Approximation of the density modulation cost. For a diffeomorphism (j) the push forward condition ) = 0£A 

can be expressed as 

0 = det {D(f))0 o <p. 


As an approximation of the source term z = d t 0 + div(ud) we take into account 

f£J0, o, 4>] = J - 5 \ det {DM ot-OYdx. 

Approximation of the dissipation cost. By Rayleigh’s paradigm BStr45l one derives models for viscous dissipation from 
elastic energies replacing elastic strains by strain rates. We proceed as in I BER15I . where a time discretization of the 
metamorphosis model was investigated and define 


[ W(Dct>)+e\D™cl>\ 2 dx. 

J D 

Here, W is a hyper elastic energy density and the higher order term |D m </>| 2 acts as a regularizing term for some small 
e > 0 and enforces the deformations to be in the space W m ’ 2 . We make the following assumptions on W (cf. also 
IBERT5) ): 

(Wl) W is non-negative and polyconvex, 

(W2) W(A) > ao|l°gdet^| — a± for ao> a i > 0, and every invertible matrix A with det A > 0, W(A) = oo for 
det A < 0, and 

(W3) W is sufficiently smooth and the following consistency assumptions with respect to the differential operator L hold 
true: IF(II) = 0, UIF(II) = 0 and \D 2 W(A)(B,B) = |(trH) 2 +/rtr((^t^) 2 ) for all B G R d ’ d . 

Due to the incorporation of this dissipation energy we finally define the space of admissible deformations over which we 
minimize in the definition of F[0, 0] as 

A = {4> £ W m,2 (D , D) : det (D<f)) > 0 a.e. in D, (j) = II on dD } , 

We assume that m >1 + 5 , which implies by Sobolev embedding that the admissible deformations are diffeomorphisms. 
Given these energy contributions we can define the compound energy 

F D [0,0, </>] = F^JO, <!>] + f£J0, 0, 0] + FlM. 

The following interpolation results justifies our choice of the time discrete path energy. 
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Theorem 3.1 (Consistency of the discrete path energy). For a convex domain D and a sufficiently smooth path of image 
intensities (0(i))te[o,i] with 9 > 0 a.e. in [0,1 \ x D and a sufficiently smooth family of velocities (f(f))tg[o,i] we consider 
interpolated images 9 k = 9(£) and motion fields v k = v(-jf). Then the resulting extended path energy 


K 


t?K,D \qK 
^5 ,7 r0 J 






“E 




fe=l 


( 21 ) 


wd/z </>£ = j^v k + II converges to the corresponding continuous path energy 


£s > 1 [ 9 , v ] : = J j 9\v\ 2 + ^z 2 +^L[v,v\dx&t. 


Proof We define the step size r = b. First, for the transport cost we easily get 


K Y [ I'/’fe - U| 2 ( 9 *-id£ = f \v£\ 2 9£_ 1 dx j f \v\ 2 9dxdt. 

k=i jD k =i Jd Jo Jd 

Following I WBRS11 1 the convergence of the dissipation cost follows from a Taylor expansion of the hyperelastic density 
function W by using the consistency assumptions: 

A' 

/ W(D<f>«) + e\D m <f>*\ 2 dx 

k=i Jd 

i K r -r 2 

= E / mU) + ^W / (U)(^f) + T ^ 2 W(U)(^f,^f) + 0(r 3 )+eT 2 |B m uf| 2 d.T 

r JD ^ 


= j D ^ ( trI}v f) 2 +^ tr 


Dv« + (Dv*F 


^ + O(t) + e\D m v k | 2 da; K ^ roo > J J L[v,v]dxdt. 


Finally, for the density modulation cost we use the Taylor expansions 

det i*f _ n + * (^) + <*>) - n + M + oM 

9foJ> I < = 6« + TVe«-v% + 0(T 2 ) 


and obtain 


k Y I \ det ( D ^kWk 0( t>k -0%-i\ 2Ax 

k= i Jd 


K 


k =1 


ID 


r\K _ r\K 

Vk Vk-1 


+ div(u fc )9 k + V9 k -v k +0 (t) 


i K—too 

da;- > 


| d t 9 + div(0u )| 2 da: df. 


to Jd 


□ 


4 Existence of time discrete geodesics 

In this section we assume that the assumptions of Section [3] are fulfilled and that 6 ,7 > 0. As before we assume that 
m > 1 + We will show that for given images 9a,9b € I = L 2 (H,K> 0 ) a time discrete geodesic exists. First we 
prove that F is well-posed in the sense that there is an optimal deformation between two images. 

Proposition 4.1 (Existence of minimizing deformations). Let 9,9 £ I. Then F D [9,9,<f>\ attains its minimum over all 
deformation (f> £ A. Moreover, (j) is a diffeomorphism and f~ l £ C 1,a (D) for a £ (0, m — 1 — |). 
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Proof. Step 1. First, we observe that F A is bounded from below, since 9 is non-negative by definition of X, W is non¬ 
negative by assumption (Wl) and the source term is non-negative anyway. Because of (W2) and <j> = II £ A there 
exists an upper bound for the energy F D on a minimizing sequence Following I BER15 I one observes that a 

subsequence, again denoted by (< iff), converges weakly in W m ’ 2 {D, D) to some </> £ W m ’ 2 {D,D) and for the limit 
deformation we get 0 _1 £ C 1,a {D). 

Step 2. We prove that <j> F r> \0. 0. c 6 ] is lower semicontinuous w.r.t. weak convergence in L 2 . It is sufficient to 

show that det {D<fP)9 o tjp — ^ det {D<p)9 o cj> in L 2 . Then the result follows from the weak lower semicontinuity of the 
X 2 -norm, the compact embedding of W m,2 {D, D) into C 1,a {D, D) for 0 < a < m — and results on the weak lower 
semicontinuity of polyconvex functionals IICia 88 l . By the assumption 9 £ L 2 (D) and by Step 1 we have a uniform 
L 2 -bound on det {D(jf) 9 o <pf so it is enough to prove that the expression converges in the sense of distributions. For 
77 £ Cf D {D) we have 


det{D(f> 3 ){x)9 o (fp {x)t]{x) dx = I 9{x)rj o {(jp) 1 (a;)da; and 


ID 


det(D(j>)(x)0 o (j)(x)rj(x)dx = / 9{x)p o 1 (®)dx. 


id 


id 


Since {fP) 


(f> 1 in C ' 1 the result follows by the dominated convergence theorem. 


Now, for a given discrete path {9q, ... ,9k ) £ X A+1 we consider Edefined in ( |2Tj ). By Proposition 


exists $ = {4> u ...,<t> K ) £ A K such that F,ff D [{9 0 ,... ,9 K ),{fii, ■ ■ ■ ,fix)] = E£ 7 [(0 o ,..., 0 K )\. Now we study 


4.1 


□ 


there 


E 


K,D 

< 5,7 


for a fixed vector of deformations. 


Proposition 4.2. Let 9 a, 9 b £ X, K > 2. Then for a fixed vector of deformations $ = {<j>i ,..., fix) £ A K there exists 
a unique discrete path {9 q, ..., 9k) £ X A+1 with 9 q = 9a and 9 k = 9b, he. 


E f' D [{e 0 ,...,9 K ),$] = inf V?' D [(0 a , 0 1 ,...,0 K - 1 ,0 b ),$\ 

(Qi,...,0k-i)£Z k 1 


Proof. First we see that the functional is bounded from above on a minimizing sequence by computing the energy of 

{9 b ,..., 9 b ) Gl 1 *- 1 : 

E ff[{9 A , 9 b ,..., 9 b , 9 b ), (0i, ■ ■ -, 4k)] < C(fi 1 ,..., 0 fe )(l + || 6 »a || 2 + IIM 2 + IIM 2 + \\0 B \\l) < 00 . 

Next we observe that for fixed {(f> x ,..., fix) £ A K the time discrete path energy is quadratically growing, i.e. 

Ef^KOAiOi . . .,0 k-1,0b), (01, • • • Ak)\ > Cl ^2 ll^lli 2 (D) ~ c 2 

i=l,...,K-l 


for constants ci, c 2 > 0 depending on {4>k)k- Therefore we can take a minimizing sequence {9\, . . . , 9 3 K _ x )j^n C X A 1 , 
which has because of the upper bound a weakly converging subsequence in L 2 with limit {9\,..., 9x~ 1 ) £ X A_1 . Now, 
the energy is strictly convex in 9k for all k = 1,..., K — 1, hence there is a unique minimizer in 1 K ~ 1 . □ 

Next, we can use these two propositions to prove existence of minimizers of the discrete path energy E f . 

Theorem 4.3 (Existence of discrete geodesics). Let 9 a, 9 b £ X, K > 2 be given. Then there exists (0\ ,..., 9k- 1 ) £ 
1 K ~ 1 s.t. 


E ^[{0 a ,0i,...,9k-i,9b)\ = 


inf E ?[(0 A ,h,...,0 K - lf 9 b )]. 

( e 1 ,...,e K - 1 )ei K ~ 1 


Proof. Taking 9° k = 9 b and 0^ = II to test the energy E A 7 we observe that the Eis bounded from above on 
a minimizing sequence (9{,..., 9 3 K _- [ )jgr-j- Take a minimizing sequence (9{,... ,9 J K _ 1 )j^ of the discrete path en¬ 
ergy E^ 7 [(0A) ■, 9 b )]. Due to Proposition |4. l| for every {9 \,.... 9 ^ K _ X ) there exists a family of optimal deforma¬ 
tions {(/){,..., (j) J K ) £ A k with Eff D [{9 A ,9{,...,9^^, 9 B ), {(/){,..., fi^)] < Eff D [(9 A , 6{,..., 9 : K _ 1 , 9 B ), for 

all T £ A h . As in the proof of Proposition^. 1 [there exists a subsequence again denoted with ( ! y i 


C l,a . By Proposition |4.2| we can assume (possible replacing (0{, ..., Of _ x ) and thereby 


for all k = 1,..., K, s.t. 

further reducing the energy) that {9 {,..., 9 J K _ 1 ) already minimizes the energy E^ AD [{9a, ■,9b ), {4,-■;<&)] in X K 1 
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Then Of. is uniformly bounded in L 2 by a constant C depending only on 0 A and On for k = 1 ..... K. This constant C is 
independent of the tpi due to the uniform bound of Xf. in IT" 1 ' 2 . Hence we can pass to a further subsequence satisfying 
(Of ,..., 0 J K _ 1 ) —=■ (0i,, Ok- i) in L 2 . To prove weak lower semicontinuity in L 2 of the functional, it is sufficient to 
pass to the limit in the identities 


det(D4> 3 k )(x)0 3 k o 4> 3 k (x)r](x) da; = / 0 3 k (x)i 7 o ( 4 > 3 k ) 1 (a:) dx 


Id 


id 


det(Dc/) k )(x)0k o 4>k(x)r](x)dx = / O k (x)r] o (f) k (x) dx , 


id 


( 22 ) 

(23) 


for an arbitrary C^-function r/, which follows from the C 1, “-convergence of (cpf ) 1 and the weak L 2 -convergence of 
0 3 k . For the demonstration of lower semicontinuity in the remaining terms we refer to analogous discussion in Proposition 

ED □ 

Finally, let us study in more detail the optimality conditions for (0 1 ,.. •, 0k-\) £ Z K ~ l in preparation of the later 
derivation of a numerical algorithm. At first we consider the simplified model without the constraint 0 k > 0 for k = 
1,..., K — 1. Since for fixed deformations the energy is strictly convex, there exists a unique minimizer. For each 
k = 1,..., K — 1 there are two terms in the energy where 0 k appears: 

F D [0 k ,0 k +i,(f)k+i} = j \(j>k+\ - U| 2 0fc + ^ det(D(f) k+ 1 )O k+1 o </> fe+1 - 0 k \' 2 dx + 7 F ^L.[0fc+i] > 

F D [ 0 fc_;L, 0 k , (j> k \ = J | (f> k - H| 2 0 fc_i + det (D(j) k )O k o (j> k - 0 fe _i [ 2 dx + 7 F^ ous [ 0 fe ] 

= Id ~ + I dk ~ ( det ( D 0 k ^ l0k -^ 0 ^l 2 det( D (f) k ) o fa 1 dx + 7 F° sco J(f) k \ . 

Hence, the Euler-Lagrange equation for 0 k is 

0 = |0 fe+ i - 1I | 2 - |(det(D0 fe+ i)0fc+i o cj) k+l - 0 k ) + ^(0 k - (det(H0fe) _ 1 0fc_i) o (j)- 1 ) det(D(j) k ) o 4 j ” 1 

for all k = 1,..., K —1 and a.e. x £ D. Now we define the discrete transport path X(x) = (X\ (x), X 2 (x ),..., Xk-i(x)) 
with Xl(x) = (f> i(x) and X k (x) = 4> k (X k -i(x)) and the vector 

0(x) = (0 1 (X 1 (x)), 0 2 (X 2 (x )),..., 0 K -i(X K -i(x))) . 


Then we can write the optimality conditions as 


det(D4> k+1 )0 k+1 o (f> k+1 + 0 k -i o cj) k 1 - | \4> k +i - H| 2 
0 k — - 


1 + det(D4> k ) o (j> 


-1 


From X k £ C 1,a we deduce that 


„ _ (det(D^fe + i) o X k )(O k+ i o X k+ i ) + 0 k _i o X k _i — ||X fc+1 — X k \ 

k k ~ 1 + det(.D<^fc) o X k _\ 


(24) 


for a.e. x £ D and for all & = 1,..., K — 1. This can be rewritten as a linear system A (x)O(x) = R(x), where A(x) is 
a tridiagonal matrix given by 

* / , det(Dcp k+1 ) o X k (x) 1 

A(x)fc fc + l — , w n , N y / \> X(x) kk ~ 1, A(x)fc k -1 — 

1 + det (Dcpk) o A fc _i(x) 

and R(x) = B(x) + T(x) with B(x), T(x) £ R x_1 given by 


1 + det(Dcj) k ) o X fe _i (x) 


B(x) 

T(x) fe 


Oa(x) det (Dc/ik) Q X K _i(x)0 B Q X K (x) \ 

1 + det(D^i)(x) ’ ’ 1 + det(D4>K-i) 0 X K - 2 (x) J 

S |X fc+1 (x)~A fc (x)| 2 

2 1+ det(Dcj) k ) o X k _i(x) 
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Now, the unique minimizer (0 1; ..., 0 k~\ ) C I K ] satisfies for a.e. x £ I) the derived linear system of equations and 
gives the only solution of this system. Thus A{x ) is invertible for a.e. x £ D and by solving the system we can recover the 
minimizer. In the constraint case 9 > 0 a.e. the minimization with respect to (0 \...., 9k-i) no longer decomposes into 
a linear system of equations with unknowns (0i(Xi(a;)), ..., 9k-i(Xk-i(x))) for a.e. x £ D. But the decomposition 
along the discrete paths (A'o(at),..., Xk(x)) is still applicable. Indeed, one observes that for a.e. x £ D the vector 
(^(X^ec)),..., 9k-i(Xk-i{x))) minimizes the quadratic functional 

Q(6i,...,6 K -i) = det(DXk~i)(x) f\X k (x ) - X fc _i(x)| 2 4-i + t| det(D</> k )(Xk-i(x))dk -4-i| 2 ^ 
k =i ' ' 

with 6*o = 9a(x) and 9 k = 9b{x) over all (0 1; ... ,9 k- i) £ M A_1 subject to the constraint 9k > 0 for all k = 
1 , ,K — 1. This is a simple quadratic optimization problem in R K 1 with inequality constraints. 


5 Spatial discretization 


With respect to the spatial discretization we follow the procedure already proposed in I BER15 I. We restrict to two 
dimensional images ( d = 2) and consider a regular quadrilateral grid on the two-dimensional image domain D = [0, l ] 2 
consisting of rectangular cells {C m } meTc with If being the associated index set. Let Vh be the space of piecewise 
bilinear continuous functions and denote by the set of nodal basis functions with / y being the index set of all 

grid nodes Xi. 

We investigate spatially discrete deformations : D —> D with <1> k £ V 2 = Vh x Vh and spatially discrete image 
maps Qk : D — > K with 0*. £ Vh- Given any finite element function U £ Vh we denote by U = ( U(x-i))i e i N the 
corresponding vector of nodal values. Now, we define a fully discrete counterpart E A 7 h of the so far solely time discrete 
path energy E defined in ( p~ 8 ] > as follows 

E A 7 j/t [(0 o ,...,0K)]= min E&? h [(0 O) ..., 0*r), (<h,. ■., *k)] 

<5>k&V 2 h 
®k\dD = H 


and obtain the resulting fully discrete approximation of the squared Riemannian distance 


nf 7j/1 [0 A ,0 s ] 2 = 


mm 

e 0 ,...,e K 

e 0 = e A ’ & R 


ei 


E 




0 / 


Here, E^^[(0 q,..., Ok), (3>i, ..., $iy)] is the discrete counterpart of E f^ D in ( |2T| ) obtained by the evaluation of all 

the integrals in E f^ D using third order Simpson quadrature with 9 quadrature points. Then, the resulting entries of the 
weighted mass matrix M/Jw, = (M/,[w, <£, . gi .^ with weight u ; and transformed via deformations $, d* 

are given by 

8 

°$)(^)(^ °^)(4). 

leic 9=0 

Here, the x l q are the quadrature points and the w l q are corresponding quadrature weights. In the case oj = 1 we write 
M ft [l,$,f] = M /l [$,4']. 

To compute a minimizer of the fully discrete energy E A 7 h we proceed as in the existence proof of time discrete 
geodesics in Section [4] and alternate the optimisation of the set of deformations for fixed image intensities and the opti¬ 
mization of the image intensities for fixed deformations. The optimization of deformations decouples in time. To calculate 
an optimal, discrete matching deformation for two consecutive images we use a conjugate gradient method for the fully 
discrete energy E f^ h . In practice we use the following hyperelastic energy W(D</)) = ^ ||D</>||^ + j (det D(f>) 2 — 
(n + a) log(detD(/)) — p, — | for det(D^) > 0 with fixed A = 10 and /.t = 1 and differing from the assumptions in 
Section [3]we skip the higher order term D m (f>\ 2 . Indeed, the associated regularization experimentally turned out not to 
be necessary, possibly due to the regularization by the spatial discretization. For a fixed vector of discrete deformations 
$ = ..., &k) the minimization of E 5,4 with respect to © = ( 0 i,..., Ok- i) leads, as in the spatially continuous 

case, to a linear system of equations. Indeed, we obtain as the discrete counterpart of f D | <j) k — A\ 2 9 k -i dec 


K 8 

£££< 

l leic q =o 


K 


4 (|$ fc - H| 2 0fc_i) 04) = 54 M h [\$ k - 1I| 2 , H, H]0 fc _!l, 


k =1 
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Figure 1: Optimal transport via translation for two different pairs of images, each of them with identical mass. Left: 
discrete geodesic between a pair of scaled bump maps with K = 4, 6 = 10 _1 , 7 = 10 -1 is shown, right: discrete 
geodesic between a square and a translated square with K = 4, S = 10” 1 , 7 = 10 -2 . 

with I = (1,..., 1) £ R Jjv and as the discrete counterpart of f D | det (D(f> k )O k 0 <t>k — 9k-i \ 2 dz 

EE EH (| det(D$fc)( 0 fe o - 0 fc _i| 2 ) (x l q ) 

k =1 l£lc Q —0 
K 

= Y, (M h [(det D<S> k ) 2 , $ fc , $ k ]Q k • Q k - 2M ft [det(D$ fc ), U]©fe • 6 fc _ 1 + M h [JL, H]©^ • 0^,) . 

k=1 

Hence, the resulting discretized part of depending on 0& is given by 

M h [|$ fc+1 - 1I| 2 , n, H]0 fc I + l(M h [(det D$ k ) 2 , $ fc , $ fc ] + M/Jll, ll])0 fe • 0 fc 

— j (M/,[det(Z3 ( hfc), H] 0 fc • 0 fc_i + M;j[det(D$fc +1 ), d) fe+1 , ll] 0 fc+1 • 0 fc ) . 

In what follows, we restrict to the non constraint case minimizing over intensities, which are not necessarily non negative. 
In fact, in our numerical experiments for 0 4 ,0« > 0 and with all deformations being initialized with the identity we did 
not observe negative density values in the vectors 0 k for k = 1 , A 1 . The implementation of a constraint, quadratic 
optimization method is work in progress. 

For the variation of with respect to the fc-th image Q k one obtains 

5 @fe E fYh = M h [|$ fc+1 - 1 I| 2 , II, II] 1 + |(M h [(det D$ k ) 2 , $ fc) $ fc ] + M h [H, U])0 fc 
(M^[det H] T 0fc_i + M/([det D$fc + i, ll]0fc + i) . 

As a consequence the necessary condition for 0 := (0i,..., Qk- 1 ) to be a minimizer of E*^ is a block tridiagonal 
system of linear equations A[<F]@ = R[$], where A[<F] is formed by (A — 1) x (K — 1) matrix blocks A k ,k' £ K /jvX/jv 
with 


Afe,fc_i = — M/jfdet 1I] T , A fcifc = M/ l [(det D$ k ) 2 , & k , <h fc ] + M/Jll, U], 

Afc^+i = —Mft[detD$fc + i,$fc + i, II] 


and R[<F] = B[<F] +T[$] consists of A — 1 vector blocks R*, = -f T/- £ K /jv with Bi = M/Jdet D&i, $ 1 , H] T 0o, 
B 2 = ... = B^_ 2 = 0, B K -1 = M^detl^^ll]©*, and T k = — fM h [|* fc+1 - H| 2 , II, U] T I for all k = 
A-l. 

The energy J2iei c Sjj=o H ~ + y I det D$ k Q k o - 0 fc _i| 2 ) (x l q ) is convex in Q k and strictly 

convex in Q k -\ . Hence, E^/^ is strictly convex in © and there is a unique minimizer 0 = 0[<b] for fixed <t>. This 
implies that A is invertible and therefore the resulting solution © coincides with the unique minimizer of E ff h - Numer¬ 
ically, the corresponding system of linear equations is solved with a conjugate gradient method with diagonal precondi¬ 
tioning. In addition, as an outer iteration of the numerical energy descent scheme we apply a cascadic approach starting 
on coarse grids and successively refining the grid. 

In what follows, we will discuss numerical results obtained by the proposed scheme. We start with two simple 
transport examples of image densities with identical mass. In Figure[T]the optimal transport geodesics connecting a bump 
map /( x) = exp ^(l — cr —2 — rco| 2 ) ^ \B a (x 0 ) with centre Xq £ D and radius a > 0 and its translate as well as 
a characteristic function of a square and its translate are considered for small 6 and 7. Indeed, the computed optimal 
transport constitutes of a translation. Next, we illustrate the role of the source term allowing for density modulation in 
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Figure 2: Discrete geodesic are computed for different viscosity parameter 7 between two bump maps placed on a square 
and periodically extended to M 2 . Top: 7 = 5 - 10 -4 , S = 10 _1 , bottom: 7 = 5 (S = 10 _1 ). The images are extracted 
from a discrete geodesic with K = 9. The different contributions to the resulting discrete path energy are: transport cost = 
0.0278363, density modulation cost = 0.00107356, dissipation cost = 0.0128489 (top row) and transport cost = 0.155598, 
density modulation cost = 0.00274782, dissipation cost = 0.00139578 (bottom row). 


case of 9a and 0/> given in Figure [3] as two bump maps of different size at different centre points and in Figure [4] as 
the characteristic functions of two rectangles of different size still for small 7. In Figure [2] we show the influence of the 
viscous dissipation. Picking up a test case from [ BBOOfl we consider image intensities on a square periodically extended 
to R 2 with a bump map once placed in the vertices of the square and once at the centre. We consider periodic boundary 
conditions both for the image intensities and for the motion field. The Wasserstein geodesic was already computed in 
HBB001 and we obtain approximately the same result for small 7. Indeed, the bump map at the vertices split up into 
four pieces, which are then transported separately into the centre. From the perspective of optimal transport this path is 
energetically preferable due to the shorter transport distance compared to a simple translation from the vertices into the 
centre. Obviously, this splitting of mass is expensive from the viscous dissipation perspective. Hence, for larger 7 we 
observe the simple translation. 

Next, we illustrate the role of the source term allowing for density modulation in case of 9a and On as in Figure[3]but 
now with different mass in the two bump maps. Still we impose periodic boundary conditions. For small values of 5 we 
observe a splitting of the bump maps in the corners with the outer one being blended out and the inner one being mainly 
transported into the middle, whereas for larger values of 6 we observe a blending process without significant transport. 
Furthermore, increasing the viscous dissipation parameter 7 leads as in Figure[2]to a translation of the whole bump, while 
the mass overhead is continuously faded-out. In Figure [4] the input images consists of characteristic functions of two 
rectangles of different size. Now, we impose natural boundary on dD. For small S and strong penalization of sources the 
surplus of mass is pushed outwards, whereas for large 5 one observes a simple blending and almost no transport. 

Furthermore, Figure [5] compares our model with the metamorphosis model on the discrete geodesic between two 
images consisting of a light and a dark square and the flipped configuration. For very small density modulation parameter 
(S = 0.01) we observe a transport of a ’’light block” from the bottom to the top square, especially mass is approximately 
preserved. In case of the metamorphism model with purely viscous flow, we see a transport of the lighter square combined 
with a fading in and out of the darker phase. 

As a first imaging application we pick up in Figure [6] an example from [ PPQ14 1. For small S and small 7 we obtain 
a very similar result. Finally, in Figure [7] the geodesic between two different slices of the same human brain recorder via 
MRI is shown. The corresponding image intensities are characterized by substantially different masses. In fact, it is the 
incorporation of both the source term and the viscous dissipation term which enables a reasonable morph between the two 
slices. Thereby, the source terms allows for local image intensity modulation, whereas the viscous dissipation ensures 
regularity of the resulting transport path. 

There is no guarantee that the alternating algorithm converges. To demonstrate the experimental convergence be¬ 
haviour we choose the application shown in Fig. [7] and show the evolution of the i 2 norm of the difference between 
consequitive intensities in Fig. [5] 


15 









Figure 3: Discrete geodesics (K = 9) between two bump maps of different mass for different values of <5 (from the first to 
the third row <5 = 1,10,100 with 7 = 5- 10 -4 ). In the fourth row the discrete geodesic for 7 = 1, <5 = 10 _1 is displayed. 



Figure 4: Discrete geodesics between two rectangles of different mass for different values of <5. Top: <5 = 10 2 , middle: 
<5 = 10“\ bottom: <5 = 1 (K = 9, 7 = 10 -2 ). 


Figure 5: Comparison of the combined model with viscosity parameter 7 = 1 (top) and the metamorphosis model 
(bottom) for <5 = 10~ 2 . 



Figure 6 : A discrete geodesic between images of Monge and Kantorovich with <5=10 2 , 7 = 10 2 (image provided by 
G. Peyre). 
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Figure 7: Two slices of the same 3D MRI data set of a human brain are connected with a discrete geodesic (data courtesy 
of H. Urbach, Neuroradiology, University Hospital Bonn). Top: discrete geodesic with S = 10 -2 , 7 = 10 _1 , bottom: 
corresponding values of Zk = det (D(j)k)0k ° (f>k — Ok -1 (blue: positive, red: negative). 



ITERATION 


ITERATION 


ITERATION 


ITERATION 


Figure 8: The convergence of the alternating descent method is shown for the application in Fig. [7] For the different levels 
of the cascadic descent scheme (K = 2,4, 8,16) the l 2 norm of the difference of consequitive space time densities (P is 
visualized using log-log plots. 
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6 The Benamou-Brenier discretization for the non viscous model 


In this section we numerically compare the proposed approach ( [TT| with the numerical scheme for optimal transport 
proposed by Benamou and Brenier l lBBOO I. where the mass constraint is relaxed. After the change of variables (9, v) i—>■ 
( 9 , m = Ov ) the minimization problem of the discrete path energy is rewritten as 

sup min ( F(B0) + G(0) + [ [ 0 ■ z — \\z\ 2 dx At 
z \ Jo Jd ° 

where 0 is a Lagrange multiplier introduced to satisfy the condition on z, F is the indicator function of the convex set 
K = {(a, b) € K. x : a + < 0}, G(0) = f D 0( 0, -)9q — 0(1, -)9\ da;, and B : 0 i-> ( d t 0 , V^). For the outer 

maximization in z one gets the optimality condition z = |</>. The augmented Lagrangian is given by 

L r [0,q,n]=F(q) + G(0) + J J 0 -z-^\z\ 2 + fi-(V ttX 0 -q)+ ^\V t , x 0-q\ 2 Axdt, 

with variables q = (a, b), jj = ( 8 , m) and Benamou and Brenier propose an alternating gradient descent to compute the 
saddle point. Using the fact that z = ^0, one updates z and 0 simultaneously solving —rA tjX 0 n + 1 0 n = di\r t x (i_i n — 
rg n_1 ) with Neumann boundary conditions in time, i.e. rd t 0 n ( 0, ■) = 6 q — 9 n ( 0, •) + ra n-1 (0, •), rd t 0 n ( 1, ■) = 9 1 — 
9 n ( 1, ■) + ra n_1 (l, •). Let us emphasize that in I B BOO I the second term on the left hand side which reflects the source 
term already appeared in the original scheme by Benamou and Brenier as a regularization term. 

To study the impact of the parameter S we pick up the problem already presented in Fig. [2] Now, we choose two input 
bump maps of different mass. Figure [9] shows discrete geodesics for different <5. For large 5 we basically observe pure 






Figure 9: The Benamou-Brenier discretization applied to optimal transport with relaxed mass constraint for input images 
with bump maps of different mass. The rows show equidistributed time steps for a time step size r F anc | different <5 
(top: (5 = 1, middle: <5 = 10, bottom: 5 = 100). 

blending and almost no transport, whereas for smaller S mass is first reduced for each bump map leading to a concentration 
in 4 bumps which are then transported. The differences to Fig. [3] seem to be due to the presence of still some viscous 
dissipation. 


7 Application of the variational time discretization to Riemannian barycentres 

As a further application of our time discrete geodesics in the space of images we consider the computation of (weighted) 
discrete barycentres. We call 9 X the barycentre of M input images 9 1 ,, 6 AI for given weights A 1 ,..., X AI with A m > 0 
and J 2 m=i ^ m = 1 > ^ minimizes 

M 

®sM9] = A m W,5, 7 [M m ] 2 

m= 1 
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Next, replacing the time continuous path energy by the time discrete energy W ’/ < 7 ( [T9| we ask for a minimizer of 
the energy 

M 

Bs-yW] = min 

( e r)f==o,...,icci 

(4> 1 k)k=l....,K<ZA 

9=e™ 

over M discrete image paths (0™)k=o,...,K (m = 1,..., M) and M discrete families (<f>™)k=i,...,K (m = 1,..., M) with 
the last image of the m discrete image paths being the ruth input image (O"} = O'" ) and the additional constraint that the 
set of first images being all equal to 6 ( 0 = 0™ for all to = 1,..., M). 

The necessary conditions for the images 0™ and the deformations oiJI' (k = 1,..., M, m = 1,..., M) are identical 
to those for simple discrete geodesics connecting the corresponding pair of images (0, 0 m ). Solely the condition for the 
barycentre image itself changes to 

M / c 

0, £ x m det- H| 2 

m—l ' 

Finally, we take into account the spatial discretization introduced in Section[5]and define 0 A as the fully discrete, weighted 
barycentre of the input images 0 1 ,..., 0 M , if 0 A minimizes the energy 




E- 

m—l 


JK r r\m 

^<5,7 ro 5 • ■ 


1 j ■ 


M 


\mi2 


= E 

m=1 

M 

min £ A m E^[0™,...,0^,$r,...,$ 


(25) 


{&k)k=0,...,KCX 

W)‘=i. 

0 = 0 ^ 


m—l 


Again for fixed deformations and skipping the non negativity constraint for the densities one obtains a 

m=l,...,M 

system of linear equations to be solved for with 0j = ... = 0 q 7 = 0 A . This linear system consists 

m— 

of M copies of the equations for (0 1; ..., &k-i) i n the system, where we replace 0^, by 0™ and <f>fe by and an 
additional set of equations for 0 A , i.e. 


M 

E 

m—l 


M h [n, H] 0 A = V A m M h [det(D$ 


l ), H]©r - 2 M h [\- n| 2 , II, 1I]I 


Still, a slightly modified strict convexity argument proves that the energy E?) h is strictly convex in the images 0™ 
for A; = 1,..., K, m = 1 ,,M and in the additional image 0 A . In particular, there exists a unique solution of the 
linear system. Let us remark that this is no longer clear if we replace /J0 m , 0] by ^[0, 0 m ] in the definition 
of the fully discrete barycenter in ( |25| ). In the implementation, we apply an analogous alternating descent scheme as 
described in Section [5] to compute fully discrete approximations of the weighted Riemannian barycenter. Furthermore, 
we use a cascadic approach, starting with coarse time discretizations and then successively refining the discretization in 
time. Figure|To|shows barycenters (with equal weights A = for three different sets of sugar beet slices extracted from 
noninvasive 3D MRI images at different days after plantation for different viscous dissipation parameters 7. Furthermore, 
we show the variability of the different contributions to the path energy between the barycenter and the input images 
for all input sugar beets. Finally, we display in Figure [IT] weighted barycenters of three different wood textures with all 
admissible combinations of A m £ {0, 1}. 
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Figure 10: Barycentres of different sets of sugar beet slices (left) for 5 = 10 -1 and for 7 = 10 —2 , 10~ 1 , 1, 10 (from left 
to right with 0j t corresponding to 7 = 10 J ). On the right the transport cost (T), density modulation cost (Z), and viscous 
dissipation cost (V) are plotted for all input slices in the case 7 = 1 (cost values for the second beet at day69 are excluded 
as outliers), (data provided by research network CROP.SENSe.net) 



Figure 11: Weighted barycentres of three different wood textures (http://de.wikipedia.org/wiki/Holz) are shown for 5 = 
10 _1 , 7 = 1 (Left: barycentric triangle where (Ao, Ai, A2) are overlaid each texture, right: discrete geodesics for K = 4 
between the three input textures and the barycentre are depicted. 
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8 Conclusion and Outlook 


In this paper we have developed a combined optimal transport and metamorphosis model and propose an effective time 
discretization of the path energy in the space of density maps. The method allows us to approximate the original Wasser- 
stein distance and for larger viscosity parameter interesting additional effects can be observed. In particular in applications 
to images the incorporated source term turns out to be an appropriate way to deal with mass variability. Let us briefly 
comment on limitations and possible future extensions of the model. So far, in the non-viscous case the source term has 
to be absolutely continuous with respect to the Lebesque measure (cf. Sectional, because the measure in the decom¬ 
position of the source measures is not unique. Therefore, a singular part in Lr would depend on the decomposition. An 
alternative model with a source term in L 1 including singular parts is work in progress. In addition, in the time discrete 
model discussed in Section[3]the treatment of the source term in L 2 required special care, since we aim at measuring the 
change of densities, which is an if -concept. Furthermore, for our generalized model including dissipation existence of 
geodesics in the time continuous case is unclear. In the non-viscous case we made use of a change of variables by consid¬ 
ering the momentum instead of the velocity, but for the viscous dissipation term this does not appear to be the appropriate 
concept. This also renders the verification of T-convergence more difficult than in the case of the metamorphosis model 
in feER15l . 
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